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1 Introduction. 

In the framework of charge transport in semiconductors, a technique widely used in or- 
der to find approximate solutions of the Boltzmann transport equation (BTE) is based 
on a spherical harmonics expansion (SHE) of the distribution function (Rahmat, White 
and Antoniadis, 1996; Vecchi and Rudan, 1998; Ventura, Gnudi and Baccarani, 1995; 
Liotta and Struchtrup, 2000). Recently an asymptotic solution of the SHE equations 
was found (Liotta and Majorana, 1999) in the case of a homogeneous (bulk) device 
with a simple parabolic band structure. Despite the very simple situation in which this 
solution was obtained, it has revealed to be very useful in order to develop new asymp- 
totic hydrodynamical models describing the hot electron population in silicon devices 
(both in the homogeneous and non-homogeneous case). In particular see Anile and 
Mascali (2000) and Anile, Liotta and Mascali (2000), where this asymptotic solution 
was used in order to close the set of moment equations. 

The aim of this work is to show the possibility of finding a new asymptotic solution gen- 
eralizing that derived in Liotta and Majorana (1999) to the case of a non-parabolic band 
structure (Kane model), this solution reduces to the old one when the non-parabolicity 
parameter goes to zero. The importance is due to the fact that the Kane equation fits 
better the real band structure in the high field regime. Therefore this solution can be 
very useful in order to develop improved high field hydrodynamical models which could 
describe better the hot electron population. This solution is also interesting by itself 
in the framework of SHE models. 

2 Basic equations. 

We consider the case of unipolar semiconductor devices in which the current is essen- 
tially due to electrons (but the results can be generalized to holes). The semiclassical 
description of the electron transport is based on the BTE (Markowich et al., 1990; 
Ferry, 1991; Cercignani, 1987), which writes 

Of t 

-± + V (k) • V x / - -E • V k / = Q(f) , (1) 

here /(£, x, k) is the electron distribution function, generally depending on time t, 
position x and wave vector k (belonging to the first Brillouin zone B). e is the absolute 
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value of the electron charge, h the reduced Planck constant, E the electric field. V x 
and Vk denote the gradient with respect to x and k respectively. The group velocity 
v(k) is determined by the conduction band structure: v(k) = ^Vke(k), where e(k) 
is the electron energy which depends on the wave vector. Q is the collision operator 
which in the non-degenerate case has the form 

Q(f) = J W(k,k)/(*,x,k)dk-/(t,x,k) J W(k,k)dk, (2) 

W(k, k) representing the electron scattering rate from a state with wave vector k to 
one with wave vector k. 

We will consider a stationary and homogeneous situation (bulk device), neglecting the 
Poisson equation and taking into account only a constant externally applied electric 
field, then it will be / = /(k). Moreover, we will suppose that the electric field 
is directed along the x-axis so to have a cylindrical symmetry around this axis and 
represent the two-dimensional momentum space by means of the polar coordinates 
k = |k| (= and 8 = arccos (k • / 1 k 1 1 | ) , where k x is the projection of k along 
x. Therefore the distribution function can be expanded in Legendre polynomials of the 
angle 9 (Rahamat et al., 1996; Liotta and Struchtrup, 2000) 

/(k) = £/ n (e)P n (co60), (3) 

n 

where P n is the nth order Legendre polynomial. This expansion will be computationally 
viable only if few spherical harmonics are enough to accurately represent the momentum 
space distribution. We will assume that the first two terms of the previous expression 
give a good approximation 

/(k)~/ o (e) + /i(e)cos0. (4) 

The lowest order harmonic coefficient furnishes information about the isotropic part 
of the distribution function and J fo dk yields the electron concentration. The first 
order harmonic coefficient describes the asymmetry of the distribution function in the 
direction of the applied electric field, and / f\ cos #v(k) dk/ J fo dk gives the hydrody- 
namical velocity of the electron gas. 

We will assume a spherically symmetric band structure of the Kane form (Ferry, 1991, 
Jacoboni and Lugli, 1989; Tomizawa, 1993) 

l{e):=e{l + ae) = (5) 

where m* is the electron effective mass and a the constant non-parabolicity parameter. 
By putting a = one obtains the usual parabolic band approximation. With this choice 
we can assume for the first Brillouin zone B = IR 3 and we have v(k) = m *(^ae+i) 1 
As regards collisions, we will take into account the interaction between electrons and 
non-polar optical phonons and that between electrons and acoustical phonons, the 
latter in the elastic approximation, valid when the thermal energy is much greater than 
that of the phonon involved in the scattering. We consider the electron scatterings 
with ionized impurities to be negligible, i.e. we assume the doping density to be low. 
Then the transition rate of the collision operator reads (Jacoboni and Lugli, 1989; 
Tomizawa, 1993) 

W(k,k) = % op [n op 5(e - e-fiuj op ) + (n op + 1) S(e -e + hoj op )) + 3C ac (5(e - e), (6) 
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where e = e(k), e = e(k), n op = (exp ( ^bTl ) ~ ~0 * s ^ e thermal equilibrium optical 
phonon number and X op and X ac are respectively the non-polar optical and acoustical 
kernel coefficients (constant at a first approximation). fiuJ op is the optical phonon 
energy, ks the Boltzmann constant and Tl the lattice temperature. These choices are 
appropriate for silicon devices. 

The SHE equations are easily obtained by inserting the expansion (||) into the BTE 
([[]) and balancing the terms of the same order in .P„(cos 9). To generate a closed set of 
equations, all coefficients of order higher than the first are set to be zero, see Rahmat 
et al. (1996) (a closure inspired by the Grad moment method, see Grad, 1958). 
But for the aims of this paper it is preferable to perform a change of variables and write 
down a set of two coupled equations in the unknowns 

N(e)=a(e)f (e), (7) 



p( £ ) = -.y 7 (£)/i(4 (8) 



where 
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a(e) := 6(e(k) - e) dk = 4^2vr ^-J H(e)( 7 (e))* 7 '(e) (9) 

is the density of states . H(e) is the Heaviside step function and ~f'(e) = = (l + 2ae). 
So doing, the expansion (||) writes 

These new variables have also a direct physical interpretation: 

/ N(e)de= / / (e)dk , / P(e)de= v(e)f 1 (e)dk, (10) 
Jo Jr 3 Jo Jr 3 

(v(e) = \v x \) and are very suitable for our problem. 

With these choices the equations of our SHE model, in the stationary homogeneous 
case, write: 

dP 

-cE— = G 1 (N) (11) 
de 

_ tE d{g{e)N) + tEh , £)N = Q , p) (12) 
de 



where 



and 



Gi{N) = X op a(e)[(n op + l)N(e + huj op ) + n op N(e-hLo op )]- 

X op [n op a(e + huj op ) - (n op + l)a(e - %u op )\ N(e) , 
Gi{P) = -[n op X op a(e + hu op ) + (n op + l)X op a(e - huj op ) + X ac a(e)]P(e) 

»M := , h(E) r 1 4 « 



3m*(7'(e)) 2 ' m* r y'{e) 3m*(y(e)) 3 

We would like to underline that equations (|iT|)-(|l2|) can be obtained directly from the 
BTE by using a new alternative procedure. It consists in multiplying both sides of 
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equation (|T|) respectively by 5(e(k) — e) and by v(k)<5(e(k) — e) and then formally 
integrating with respect to k over the whole space M 3 . Some suitable closure relations 
are needed: in particular by assuming that / depends on k only through the variable 
e one obtains equations (|i"l|)-(|T2|) in the general non-homogeneous, non-stationary case 
(Liotta and Majorana, 1999; Majorana, 1998). This method is similar to the method 
of frequency dependent moments of radiation hydrodynamics (Thorne, 1981). 



3 Dimensionless equations and physical parameters. 

It is useful to introduce dimensionless variables: let 



is* 



4\/2- 



7T 



\/k B T L X 



op^op 



V m* J 



w:=—, n(w) :=uJlN(e), p(w) := u*{®UP(e) , 



* fo^op 



a := — = e , k :- 

Hop 



, C : = lE— , (3 := a e. 



flop X-op ^* 



Moreover in the following we put x( w ) '■= + P w ) an d x'( w ) = = (1 + 2j3w). By 



using these new variables, equations (0)-(H) become 
dp 

(- — = u(w) [an(w + A) + n(w — A)l — [u(w + A) + au(w — A)l n(w) 
dw 



d(r(w)ri) 
dw 



where 



+ q(w)n 

fj,(w) 
r(w) 

q(w) 



[fj,(w + A) + aa(w — A) + Kfj,(w)] p{w) 



H{w) [xH) 2 x'H 

2 xH 



3 [x'M] 
1 



4 X(w) 



(13) 
(14) 



X'(w) ^[x'(w)} 3 ' 
We associate the following conditions to equations (|T^)-([l^) 
n(0) = (=► p(0) = 0) , lim p(w) = 

w— >+oo 
r+00 

n{w) > Vw>0, / n(w) dw > and < +00. 

J 

Since equations (|l^)-([l^) are linear and homogeneous, if a solution (n(w),p(w)), satis- 
fying the above conditions exists, then, for every c > 0, also (cn(w) , cp(w)) is a solution. 
The appropriate values of the physical parameters, in the case of a silicon device, are 
given in table I, were m e denotes the electron rest mass. 
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7TJ * — n m 


Tt — 


300 k 


T)ii> — 063 pV 

/ VLU Q'p U . \J\J\j C V 


op 87r 2 pa; p 


D t K 


o 

= 11.4 eV A -1 


p = 2330 Kg m~ 3 




= 


9 eV 


vq = 9040 m sec -1 . 


a = 0.5 eV- 1 









Table I. Values of the physical parameters used in this paper. 
Using these parameters, we get A ~ 2.437, n ~ 5.986 and ~ 0.012926. 



4 Asymptotic equations. 

Now we want to show that it is possible to find an approximate solution of the equations 
([h|)-(|l4|) valid for high values of the electron energy. 
It is useful to introduce a new variable ip(w) defined by 



Equations (|!3|)-(|l4|) become 
dp 



-c 



dw 



3 x'( w ) dw 



n(w) = n(w)ip(w) 



n(w) [a/i(w + \)tp(w + A) + fi(w — X)ip(w — A)] 
n(w) [fi(w + A) + afi(w — A)] tp(w) 

— [n(w + A) + afi(w — A) + Kp(w)] p(w) . 



(15) 



(16) 
(17) 



Because we look for an asymptotic form of the equations (|T^)-(|T^) for large values of 
the energy w, we expand the coefficients of the equations up to the zeroth order in A: 
ji{w ± A) ~ p{w). In this way we obtain a new set of equations 



c 



dpA 



dw 
Pa{w) 



H 2 {w) [aipA(w + A) + i/ja(w - A) - (a + 1)^a(w)] 



[xH1 ; 



dtp a 



1(- 

3 x'{ w )^{ w ) [1 + a + t] dw 



(18) 
(19) 



where the subscript A label the new unknowns. Substituting (19) into (jig), it follows 



c 2 



3(1 + o + k) 



X(w) „ 



1 4faM 

x'H [x'HY 



X(w) [x(w)] 2 (aifjA (w + A) + ip A (w - A) - (a + l)ip A (w)) 



(20) 



where the primes denote derivatives with respect to w. 
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5 Approximate solution. 



In order to find an approximate solution of equation (^) we expand the coefficients up 
to the first order in the non-parabolicity parameter (3. 



w - 3f3w 2 ) $' A + (1 - 6pw) if/ A 



2 C 2 

3 (1 +CL + K) 

(w + 5f3w 2 ^) (aip A (w + A) + i\> A (w - A) - (a + l)ij> A (w)) 



(21) 



This choice is justified by the smallness of (3 and by the Kane equation itself, which is 

of the first order in the non-parabolicity parameter. 

We will search for solutions of equation ( |2~l| ) having the form 

(22) 



iPa(w) = e^ w ^ , with f{w) := t]qw + r]i(3w 2 



where r/o and n\ are functions of the applied electric force £. 

It is useful to observe that f(w±\) = f(w) + f(±\)±2r}iP\w. Expanding the following 
quantities up to the first order in (3 : 



±2 m /3Xw 
P /(±A) 



1 ± 2r)xP\w, 



substituting (^) into ( |2l| ) and retaining only terms up to first order in (3, we obtain, 
after dropping the common factor e^ w \ the equation 



c 2 



Vo + (vo + 4 /fyi ~ 6/3%) w + (-3Pr$ + 4/3r/ r/i 



3(l + a + k) 

(a e* A + e- r '° x ) (l + ^A 2 ) - (a + 1) 



tu- 



rn + 



a e noX (2rn/3\ + 5/3) + e~ voA (-2rn0\ + 5/3) - 5/3(a + 1) 



-?7oA 



w 



(23) 



If we divide both sides of equation ( p3|) by w 2 , and neglect the terms in but not 
those in — (in some sense we are serching for a " weakly asymptotic" solution), we obtain 
the following system of two transcendent equations in the unknowns 770 and rji 



c 2 



3(l + a + 
2 C 2 



— (if - 6/?77o + 4/%) = (a e* A + e~* A ) (l + ^A 2 ) 



fa + 1) 



(24) 



3(1 + o + k) 



4r7 r/i - 37$ ) = a e wA (5 + 2r/iA) + e"' 70 '* (5 - 2t/iA) - 5(o + 1) , (25) 



10A 



where in the second equation we have dropped the common factor (3. If we are able to 
solve the previous system, it is possible to obtain 770 and 771 as functions of the applied 
electric force £, and then tp A {w). Moreover, by using equation (IE), one can find pa{w)- 
Henceforth we will indicate as asymptotic solution of the SHE equations the expressions 
of ua anp pa which can be obtained by means of the approximate solution of (|l8|)-(pl 
which have been found. 
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6 Discussion of the solution and comparison with numer- 
ical results. 



An analytical solution of the equations (g3)-(ga) has turned out to be impossible. 
Therefore we have limited ourselves to a graphical and numerical analysis. Given a 
value of the applied electric force £, the requirement that the functions be integrable in 
[0, -|-oo[ tells us that both j]q and r}\ have to assume negative values. In fact, we found 
a negative solution of the system (|24|)-(|25|) in the domain [—1,0] x [—1,0] of the (r/o, 
plane, for all the values of the electric field in the explored range. We used a simple 
MapleV algorithm in order to find the solutions. In table II we give some of the values 
we found. 



E (V/cm) 


Vo 


m 


0.0 


-1.0 


0.0 


1.0 x 10 2 


-0.9999636274 


-0.0001493648 


1.0 x 10 3 


-0.9963693066 


-0.0148712049 


5.0 x 10 3 


-0.9136708716 


-0.3294754506 


1.0 x 10 4 


-0.7162321385 


-0.8333494338 


3.0 x 10 4 


-0.2511364699 


-0.5681222515 


5.0 x 10 4 


-0.1270766899 


-0.2772663271 


7.0 x 10 4 


-0.0780599266 


-0.1591756496 


1.0 x 10 5 


-0.0452815087 


-0.0850850173 



Table II. Some values of r]o and rji as functions of the electric field. 

In figures 1, 2 and 3, we compare the asymptotic solution (jia,Pa) and the numerical 
solution (un,Pn) of equations (H)-(0), the latter obtained by a suitable numerical 
technique (Liotta and Majorana, 1999), for some values of the applied electric field. 
As is possible to see, the agreement between the numerical and asymptotic solutions 
is good in all the energy range, despite the obtained asymptotic solution should be 
good only for high enough energy values (the agreement is obviously very good in this 
region). We want to observe that for low values of the electric field (|E| < 10 3 V/cm), 
instead of formula (jlfi|), we use the following expression for pa- 



1 



PA(w) 



fj,(w + A) + an(w 



2 [xHIi^m 

A) + Kji{w) 3 x'( w ) dw 



(26) 



Electric Field = 10 3 V/cm Eleclric Field = 10 3 V/cm 




Energy, (phonon energy = 1| Energy, (phonon energy = 1) 

Figure 1: Electric field = 10 3 V/cm. n and p versus energy. Dots: numerical solution, 
continuous line: asymptotic solution. The optical phonon energy is set equal to one. 
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This choice is not coherent with the expansion, but allows a better agreement with the 
numerical solution and puts in evidence a discontinuity in the derivatives at the point 
w = A (in dimensional variables e = huj op ) due to the term fj,(uu — A), that is zero for 



< w < A. This term is also present in the original set of equations (|13|)-(14). Then, 
such discontinuity is also expected in the solutions. 

As a measure of the difference between the numerical and asymptotic solution we 
can also compute V as = J pAde/ J nAde and V num = J pNde/ J riNde, which give, in 
dimensional units, the electron hydrodynamical velocity. In table III we compare the 
values we found. 



E (V/cm) 


V as (m/sec) 


Vnum(m/sec) 


1.0 x 10* 


1.4879 x 10 3 


1.4537 x 10 3 


1.0 x 10 3 


1.4831 x 10 4 


1.3858 x 10 4 


5.0 x 10 3 


7.0223 x 10 4 


5.1060 x 10 4 


1.0 x 10 4 


6.0652 x 10 4 


7.3611 x 10 4 


3.0 x 10 4 


9.8325 x 10 4 


9.7714 x 10 4 


5.0 x 10 4 


1.0106 x 10 5 


9.9872 x 10 4 


7.0 x 10 4 


9.5124 x 10 4 


9.8395 x 10 4 


1.0 x 10 5 


8.2528 x 10 4 


9.5101 x 10 4 



Table III. Comparison between the electron hydrodynamical velocities calculated by using 
respectively the asymptotic and numerical solutions. 



The behaviour of the hydrodynamical velocity is the typical one when the Kane equa- 
tion is used (see Tomizawa, 1993, p. 100, fig. 3.11). We do not consider higher values 
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of the electric field because in this case also the Kane equation becomes inadequate to 
describe the conduction band, and a full band structure should be used (Vecchi and 
Rudan, 1998). 

7 Conclusions and acknowledgments. 

We have shown that it is possible to find an " asymptotic solution" of the SHE equation 
(^)-(|IJ) in which the dependence on the applied electric field is given implicitly through 
the system of transcendent equations (^4|) -(^5|). The solution of the system (|24[)-(|25l) 
can be obtained by using standard numerical techniques. The agreement with the 
numerical solution is good in all the explored range of applied electric fields. 
If we put [5 = 0, we recover the parabolic band case and the asymptotic solution 
reduces to the asymptotic one already found by Liotta and Majorana (1999) (L-M 
solution). The L-M solution was used by Anile and Mascali (2000) to obtain a two 
fluid hydrodynamical model, where the electron population is splitted in two sub- 
populations: low-energy and high-energy electrons, separated by a threshold energy. 
In order to close the moment equations relatively to the hot electrons they used a 
distribution function whose form is directly suggested by the L-M solution, whereas 
relatively to cold electrons a a standard maximum-entropy distribution function is 
utilized. A detailed study of the resulting high-energy hydrodynamical model is given 
in Anile, Liotta and Mascali (2000). Some works about the extension of this model 
to the case of Kane band, using the asymptotic solution found in this paper, are in 
progress. A better description of the hot electron population is expected. 
The author acknowledges support from Italian CNR (Prog. N. 96. 03855. CT01), from 
TMR (Progr. n. ERBFMRXCT970157), from Italian MURST (Prot. n. 9801169828- 
005) and from " Convenzione quadro Universita di Catania-ST Microelectronics" . 
The author wish also to thank prof. A. M. Anile, prof. A. Majorana and dr. G. Mascali 
for useful discussions and precious suggestions. 
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